Getting started with Sentinel 2 data

Background As of mid-February 2018, Sentinel 2 data is available to access within the a development environment on AWS. There are a number of things that need to be done prior to gaining access to the Sentinel 2 archive. For the purpose of this notebook, we will assume you have successfully gained access to the AWS environment where Sentinel 2 data is currently housed.

What does this notebook do? This notebook steps you through how to load in and plot up data from Sentinel 2. It explores the data that are available from the Sentinel 2 satellite, and briefly describes the Sentinel satellite bands. It then loads in the s2a_ard_granule product and plots it up in true and false colour. It uses the provided pixel quality data to filters the example scene based on clear pixels.

Date: February 2018.

Author: Claire Krause

Tags Sentinel2, products, threeBandImage, threeBandImage_subplots, dc.list_products, dc.load, query, beginner, plot, image, pixelquality

% pylab notebook
from datacube.storage import masking
from datacube import Datacube
from datetime import datetime
from skimage import exposure

# Replace '156' and 'cek156' with the path to your own home directory on the VDI
s2aws = Datacube(config='/home/156/cek156/.aws_datacube.conf')

Populating the interactive namespace from numpy and matplotlib

Define some plotting functions for later on

def threeBandImage(ds, bands, time = 0, figsize = [10,10], projection = 'projected'):
    '''
    threeBandImage takes three spectral bands and plots them on the RGB bands of an image.

    Inputs:
    ds -   Dataset containing the bands to be plotted
    bands - list of three bands to be plotted

    Optional:
    time - Index value of the time dimension of ds to be plotted
    figsize - dimensions for the output figure
    projection - options are 'projected' or 'geographic'. To determine if the image is in degrees or northings
    '''
    t, y, x = ds[bands[0]].shape
    rawimg = np.zeros((y,x,3), dtype = np.float32)
    for i, colour in enumerate(bands):
        rawimg[:,:,i] = ds[colour][time].values
    rawimg[rawimg == -999] = np.nan
    img_toshow = exposure.equalize_hist(rawimg, mask = np.isfinite(rawimg))
    fig = plt.figure(figsize = figsize)
    imshow(img_toshow)
    ax = plt.gca()
    ax.set_title(str(ds.time[time].values), fontweight = 'bold', fontsize = 16)
    ax.set_xticklabels(ds.x.values)
    ax.set_yticklabels(ds.y.values)
    if projection == 'geographic':
        ax.set_xlabel('Longitude', fontweight = 'bold')
        ax.set_ylabel('Latitude', fontweight = 'bold')
    else:
        ax.set_xlabel('Eastings', fontweight = 'bold')
        ax.set_ylabel('Northings', fontweight = 'bold')

def threeBandImage_subplots(ds, bands, num_cols, figsize = [10,10], projection = 'projected', left  = 0.125,
                            right = 0.9, bottom = 0.1, top = 0.9, wspace = 0.2, hspace = 0.4):
    '''
    threeBandImage_subplots takes three spectral bands and multiple time steps, and plots them
    on the RGB bands of an image.

    Inputs:
    ds -   Dataset containing the bands to be plotted
    bands - list of three bands to be plotted
    num_cols - number of columns for the subplot

    Optional:
    figsize - dimensions for the output figure
    projection - options are 'projected' or 'geographic'. To determine if the image is in degrees or northings
    left  = 0.125  # the space on the left side of the subplots of the figure
    right = 0.9    # the space on the right side of the subplots of the figure
    bottom = 0.1   # the space on the bottom of the subplots of the figure
    top = 0.9      # the space on the top of the subplots of the figure
    wspace = 0.2   # the amount of width reserved for blank space between subplots
    hspace = 0.2   # the amount of height reserved for white space between subplots
    '''
    # Find the number of rows/columns we need, based on the number of time steps in ds
    timesteps = ds.time.size
    num_rows = int(ceil(timesteps/num_cols))
    fig, axes = plt.subplots(num_rows, num_cols, figsize = figsize)
    fig.subplots_adjust(left  = left, right = right, bottom = bottom, top = top, wspace = wspace, hspace = hspace)
    numbers = 0
    try:
        for ax in axes.flat:
            t, y, x = ds[bands[0]].shape
            rawimg = np.zeros((y,x,3), dtype = np.float32)
            for i, colour in enumerate(bands):
                rawimg[:,:,i] = ds[colour][numbers].values
            rawimg[rawimg == -999] = np.nan
            img_toshow = exposure.equalize_hist(rawimg, mask = np.isfinite(rawimg))
            ax.imshow(img_toshow)
            ax.set_title(str(ds.time[numbers].values), fontweight = 'bold', fontsize = 12)
            ax.set_xticklabels(ds.x.values, fontsize = 8, rotation = 20)
            ax.set_yticklabels(ds.y.values, fontsize = 8)
            if projection == 'geographic':
                ax.set_xlabel('Longitude', fontweight = 'bold', fontsize = 10)
                ax.set_ylabel('Latitude', fontweight = 'bold', fontsize = 10)
            else:
                ax.set_xlabel('Eastings', fontweight = 'bold', fontsize = 10)
                ax.set_ylabel('Northings', fontweight = 'bold', fontsize = 10)
            numbers = numbers + 1
    except IndexError:
        # This error will pop up if there are not enough scenes to fill the number of rows x columns, so we can
        # safely ignore it
        pass

See what Sentinel 2 products are currently available

products = s2aws.list_products()
display_columns = ['name', 'description', 'instrument', 'platform', 'product_type', 'crs', 'resolution']
sentinel_products = products[products['instrument'] == 'MSI'][display_columns]
sentinel_products

	name	description	instrument	platform	product_type	crs	resolution
id
79	s2a_ard_granule	Sentinel-2A MSI ARD - NBAR NBART and Pixel Qua...	MSI	Sentinel-2A	S2MSIARD	NaN	NaN
69	s2a_level1c_granule	Sentinel-2A Level1C - Ortho Rectified Top of A...	MSI	Sentinel-2A	S2MSI1C	NaN	NaN
71	s2a_sen2cor_granule	Sentinel-2 Level 2 - Sen2Cor Bottom of Atmosph...	MSI	SENTINEL_2A	S2MSI2Ap	NaN	NaN
81	s2b_ard_granule	Sentinel-2B MSI ARD - NBAR NBART and Pixel Qua...	MSI	Sentinel-2B	S2MSIARD	NaN	NaN
70	s2b_level1c_granule	Sentinel-2B Level1C - Ortho Rectified Top of A...	MSI	Sentinel-2B	S2MSI1C	NaN	NaN

There are two spectral Sentinel satellites currently in DEA; Sentinel-2A and Sentinel-2B. Sentinel-2A was launched on 23 June 2015 and Sentinel-2B followed on 7 March 2017.

Both of the Sentinel 2 satellites carries an innovative wide swath high-resolution multispectral imager with 13 spectral bands. The mission is based on a constellation of two identical satellites in the same orbit, 180° apart for optimal coverage and data delivery. Together they cover all Earth’s land surfaces, large islands, inland and coastal waters every five days at the equator.

For more information on the Sentinel 2 platforms and applications, check out the European Space Agency website.

Now we want to actually load and look at some data

We will focus on data from the Sentinel-2a platform for this demonstration, as there is more data available than 2b. To explore the spectral datasets from Sentinel-2a, we will use the s2a_ard_granule product. Sometimes multiple scenes are acquired by the satellite on the same day. We want to group these together, and will use group_by='solar_day' to do this.

To load in the Sentinel data, we have a number of options we can use in building our data extraction query:

• lat/lon - specify a bounding box for the extraction. Note that a polygon can be used instead. See here for more details
• output_crs - the output coordinate reference system to project the data into. The CRS you request will impact on the format for the resolution query. Two handy CRSs are ‘EPSG:3577’, which is the Australian Albers projected coordinate system, and ‘EPSG:4325’, which is WGS84 a global geographic coordinate system.
• resolution - the requested output resolution for the data. If you have selected a geographic coordinate system for the output_crs, this will be in degrees. If you have selected a projected coordinate system, this will be in metres.
• time - the time range for the query.

query = {
        'lat': (-35.25, -35.35),
        'lon': (149.05, 149.17),
        'output_crs': 'EPSG:3577',
        'resolution': (-10, 10),
        'time':('2017-01-01', '2017-02-15')
        }

Canberra = s2aws.load(product='s2a_ard_granule', group_by='solar_day', **query)

See what came back from the extraction

Canberra

<xarray.Dataset>
Dimensions:        (time: 4, x: 1229, y: 1248)
Coordinates:
  * time           (time) datetime64[ns] 2017-01-05T00:02:12.026000 ...
  * y              (y) float64 -3.953e+06 -3.953e+06 -3.953e+06 -3.953e+06 ...
  * x              (x) float64 1.542e+06 1.542e+06 1.542e+06 1.542e+06 ...
Data variables:
    pixel_quality  (time, y, x) uint8 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
    contiguity     (time, y, x) uint8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
    aerosol        (time, y, x) int16 886 886 886 886 886 931 931 931 931 ...
    blue           (time, y, x) int16 777 703 714 732 814 782 752 773 690 ...
    green          (time, y, x) int16 904 839 843 865 910 877 837 816 739 ...
    red            (time, y, x) int16 987 843 858 828 963 907 822 892 795 ...
    rededge1       (time, y, x) int16 1291 1303 1303 1282 1282 1221 1221 ...
    rededge2       (time, y, x) int16 1975 1931 1931 1945 1945 1887 1887 ...
    rededge3       (time, y, x) int16 2336 2346 2346 2337 2337 2316 2316 ...
    nir            (time, y, x) int16 2177 2628 2573 2493 2357 2135 2351 ...
    rededge4       (time, y, x) int16 2647 2637 2637 2573 2573 2537 2537 ...
    swir1          (time, y, x) int16 2223 2114 2114 2039 2039 1972 1972 ...
    swir2          (time, y, x) int16 1527 1464 1464 1410 1410 1316 1316 ...
    t_contiguity   (time, y, x) uint8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
    t_aerosol      (time, y, x) int16 879 879 879 879 879 923 923 923 923 ...
    t_blue         (time, y, x) int16 771 694 703 718 796 764 735 755 677 ...
    t_green        (time, y, x) int16 897 828 829 848 890 857 817 797 725 ...
    t_red          (time, y, x) int16 980 831 843 811 940 884 802 870 779 ...
    t_rededge1     (time, y, x) int16 1284 1281 1281 1254 1254 1192 1192 ...
    t_rededge2     (time, y, x) int16 1965 1900 1900 1903 1903 1842 1842 ...
    t_rededge3     (time, y, x) int16 2325 2314 2314 2296 2296 2272 2272 ...
    t_nir          (time, y, x) int16 2163 2598 2538 2452 2312 2094 2305 ...
    t_rededge4     (time, y, x) int16 2635 2602 2602 2528 2528 2489 2489 ...
    t_swir1        (time, y, x) int16 2211 2078 2078 1993 1993 1924 1924 ...
    t_swir2        (time, y, x) int16 1518 1438 1438 1377 1377 1283 1283 ...
Attributes:
    crs:      EPSG:3577

About Sentinel 2 bands

Sentinel 2 satellites have 13 spectral channels:

Sentinel 2 bands	DEA band name	Band number	Central wavelength (nm)	Resolution (m)	Bandwidth (nm)
Coastal aerosol	aerosol	1	443	60	20
Blue	blue	2	490	10	65
Green	green	3	560	10	35
Red	red	4	665	10	30
Vegetation red edge	rededge1	5	705	20	15
Vegetation red edge	rededge2	6	740	20	15
Vegetation red edge	rededge3	7	783	20	20
NIR	nir	8	842	10	115
Narrow NIR	rededge4	8A	865	20	20
Water vapour	N/A	9	945	60	20
SWIR - Cirrus	N/A	10	1375	60	20
SWIR	swir1	11	1610	20	90
SWIR	swir2	12	2190	20	180

These bands cover the visible, near-infrared and short-wave infrared wave lengths

Sentinel 2 bands

Data corrections

There are two corrections applied to the Sentinel data:

1. NBAR (e.g. ``green``)

   NBAR stands for Nadir-corrected BRDF Adjusted Reflectance, where

   BRDF stands for Bidirectional reflectance distribution function The approach involves atmospheric correction to compute surface-leaving radiance, and bi-directional reflectance modelling to remove the effects of topography and angular variation in reflectance.

2. T_NBAR (e.g. ``t_green``)

   Surface reflectance T_NBAR includes the terrain illumination

   reflectance correction and has the same features of NBAR, along with some additional features.

Note that the t_nbar data insert a missing value (-9999) into the dataset to denote a terrain shadow. This should be converted to a NaN value before using this data to avoid treating it as a valid value.

For more information on these corrections, see the explanation of the Landsat archive data.

Plot this scene up as true colour image

True colour images approximate what the human eye sees when looking at a landscape. Note that the function used for this true colour plot enhances the contrast between the bands, resulting in a colour-enhanced image.

threeBandImage(Canberra, bands = ['red', 'green', 'blue'], time = 1)

Plot this scene up as false colour image

This plot uses the SWIR, NIR and green bands to accentuate the presence of water in the landscape.

threeBandImage(Canberra, bands = ['swir1','nir','green'], time = 1)

Compare some scenes side-by-side

threeBandImage_subplots(Canberra, bands = ['swir1','nir','green'], num_cols = 2, figsize = [10, 10], wspace = 0.35)

Filter scenes using pixel quality

A pixel quality dataset is provided alongside the spectral data, which can be used to filter out noisy pixels. The key things we want to filter are clouds and shadows.

The pixel quality field contains values from 0 to 5.

Value	Description
0	Null
1	Valid
2	Cloud
3	Cloud shadow
4	Snow
5	Water

In the Sentinel data, clear pixels are denoted as 1 in the pixel quality mask. We can use this to filter a cloudy image.

Here is the cloudy image we would like to mask

threeBandImage(Canberra, bands = ['swir1','nir','green'], time = 3)

Plot up the pixel quality information for the same scene

Canberra.isel(time = 3).pixel_quality.plot(figsize = [10, 10])

/g/data/v10/public/modules/agdc-py3-env/20171214/envs/agdc/lib/python3.6/site-packages/xarray/plot/utils.py:51: FutureWarning: 'pandas.tseries.converter.register' has been moved and renamed to 'pandas.plotting.register_matplotlib_converters'.
  converter.register()

<matplotlib.collections.QuadMesh at 0x7fdff8340eb8>

Now use the pixel quality information to create a mask, and apply it to the spectral data

clear_pixels = Canberra.pixel_quality == 1
Clear_Canberra = Canberra.where(clear_pixels)

threeBandImage(Clear_Canberra, bands = ['swir1','nir','green'], time = 3)